N. harperi and small mammal communities
Elliott Schmidt
2025-03-26
Scenario
load libraries
library(tidyverse)
library(readxl)
library(grid)
library(lubridate)
library(janitor)
library(brms)
library(rstan)
library(reshape2)
library(modelr)
library(emmeans)
library(bayestestR)
library(tidybayes)
library(patchwork)
library(ggpubr)
library(cowplot)
library(HDInterval)
library(ggridges)
library(DHARMa)
library(rphylopic)
library(ggeffects)Data import
orangemitedata_2020 <- read_delim("data/orangemitedata_2020.txt",
"\t",
escape_double = FALSE,
col_types = cols(Date = col_datetime(format = "%d/%m/%Y"),
Trapline = col_factor(levels = c("1","2", "3", "4","5", "7", "8","9")),
Species = col_factor(levels = c("1","2", "3")),
Sex = col_factor(levels = c("1","2")),
Year = col_factor(levels = c("1","2", "3", "4"))),
trim_ws = TRUE)|>
rename(OrangeMites =`Orange Mites (Y/N)`)
smallmammaldata3 <- read_delim("data/smallmammaldata_v3.txt",
"\t",
escape_double = FALSE,
col_types = cols(Date = col_date(format = "%Y-%m-%d"),
Trapline = col_factor(levels = c("1","2", "3", "4", "5", "6", "7","8","9")),
Species = col_factor(levels = c("1","2", "3")),
Sex = col_factor(levels = c("1","2")),
Year = col_factor(levels = c("1","2", "3")),
X14 = col_skip(),
X15 = col_skip()),
trim_ws = TRUE)|>
rename(OrangeMites =`Orange Mites (Y/N)`)
abund <- read_excel("data/smallmammaldata2016-2018_all.xlsx")
abund_2020 <- read_delim("data/smallmammaldata2020_all.txt",
"\t", escape_double = FALSE, col_types = cols(Date = col_date(format = "%d/%m/%Y")),
trim_ws = TRUE) Data manipulation
orangemitedata_v2 <- bind_rows(smallmammaldata3, orangemitedata_2020) |>
mutate(Trapline =factor(Trapline, levels =c('4','1','2','3','5','6','7','8','9')),
Species =factor(Species, levels =c('2','1','3')),
OrangeMites =factor(OrangeMites),
Year =factor(Year),
TagLeft =case_when(TagLeft == NA ~ TagRight,
TRUE ~ TagLeft),
Date.Julian =as.numeric(format(as.POSIXct(Date, formate="%Y-%m-%d"), "%j"))) |>
select(c(Date:Year, Date.Julian)) |>
clean_names(case="all_caps")
abund_1 <- rbind(abund,abund_2020) |>
clean_names() |>
mutate(orange_mites_y_n = toupper(orange_mites_y_n),
species = toupper(species),
orange_mites_y_n = replace_na(orange_mites_y_n, "N"),
year = year(date),
month = month(date),
tag_left = coalesce(tag_left,tag_right),
tag_left = ifelse(tag_left=="RIP", tag_right, tag_left)) |>
group_by(year) |>
mutate(first_om_appearence = min(date[orange_mites_y_n=="Y"])) |>
ungroup() |>
subset(date > first_om_appearence) |>
mutate(month_name = month.abb[month],
species = case_when((species == "SM") ~ "DM",
(species == "NFLYER") ~ "FLYER",
TRUE ~ species),
trapline_group = as.numeric(trapline), #make a separate group that groups traplines based on habitat type
trapline_group = case_when((trapline >100 & trapline<200) ~ 100,
(trapline >200 & trapline<300) ~ 200,
(trapline >300 & trapline<400) ~ 300,
(trapline >400 & trapline<500) ~ 400,
(trapline >500 & trapline<600) ~ 500,
(trapline >700 & trapline<800) ~ 700,
(trapline >800 & trapline<900) ~ 800,
TRUE ~ 900)) |> #place new column next to the 'trapline' column to make sure it worked
relocate(trapline_group, .after=trapline) |>
unite("group_ID", year, month_name, trapline_group, remove = F)
abund_2 <- abund_1 |>
select(c("group_ID":"species", "tag_left", "year":"month_name")) |>
group_by(year) |>
distinct(tag_left, .keep_all = T) |>
ungroup() |>
reshape2::dcast(group_ID~species, length) |>
select(-c(`TRAP MALFUNCTION`, SOREX)) |> #`TRAP DISTURBANCE`
mutate(total =rowSums(across(-group_ID)),
perc_rbv = RBV/total) ## Using month_name as value column: use value.var to override.#join 'abund_2' dataframe to 'orangemitedata_v3' dataframe by the group_ID column
#so that the percentage of voles within a community (per month) can be recorded for each observation
orangemitedata_v3 <- orangemitedata_v2 |>
mutate(year = year(DATE),
month = month(DATE),
month_name = month.abb[month],
trapline_group= as.numeric(as.character(TRAPLINE))*100) |>
unite("group_ID", year,month_name,trapline_group, remove=F) |>
left_join(abund_2[,c(1,11)], by="group_ID") |>
clean_names(case="all_caps")
orangemitedata_y1y2 <- orangemitedata_v3|>
select(c(TRAPLINE, SPECIES, TAG_LEFT,ORANGE_MITES,YEAR,DATE_JULIAN)) |>
filter(TRAPLINE != 9) |>
filter(YEAR ==1 | YEAR ==2) |>
mutate(TRAPLINE =factor(TRAPLINE, levels=c('4','1','2','3','5','7','8')),
YEAR =factor(YEAR, levels =c('1','2')))
orangemitedata_y3y4 <- orangemitedata_v3|>
select(c(TRAPLINE, SPECIES, TAG_LEFT,ORANGE_MITES,YEAR,DATE_JULIAN)) |>
filter(TRAPLINE != 9) |>
filter(YEAR ==3 | YEAR ==4) |>
mutate(TRAPLINE =factor(TRAPLINE, levels=c('4','1','2','3','5','7','8')),
YEAR =factor(YEAR, levels =c('3','4'))) Descriptive statistics
temp1 <- separate(abund_2, group_ID, c("year","month","trapline.no"))|>
filter(year=="2016"|year=="2017")
summarise_all(temp1[4:12],funs(sum))/sum(temp1[12]) ## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.temp1 |> select("trapline.no","DM","RBV","WJM") |>
mutate(trapline.no =factor(trapline.no)) |>
group_by(trapline.no) |>
summarise_each(list(sum))## Warning: `summarise_each()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.temp2 <- separate(abund_2, group_ID, c("year","month","trapline.no"))|>
filter(year=="2018"|year=="2020")
summarise_all(temp2[4:12],funs(sum))/sum(temp2[12])## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.temp2 |> select("trapline.no","DM","RBV","WJM") |>
mutate(trapline.no =factor(trapline.no)) |>
group_by(trapline.no) |>
summarise_each(list(sum))## Warning: `summarise_each()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.Data analysis
High-RBV year
Fit models
priors
model formula
Fit model
model validations (convergence)
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model investigation
summary
## Family: bernoulli
## Links: mu = logit
## Formula: ORANGE_MITES ~ SPECIES + TRAPLINE + YEAR + scale(DATE_JULIAN, center = TRUE, scale = TRUE) + SPECIES:TRAPLINE + SPECIES:scale(DATE_JULIAN, center = TRUE, scale = TRUE) + (1 | TAG_LEFT)
## Data: orangemitedata_y1y2 (Number of observations: 1446)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~TAG_LEFT (Number of levels: 603)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.74 0.30 2.19 3.35 1.00 1038 1600
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI
## Intercept 7.04 1.04 5.16
## SPECIES1 -4.87 1.16 -7.22
## SPECIES3 -4.15 1.33 -6.75
## TRAPLINE1 -6.39 1.80 -9.95
## TRAPLINE2 -2.21 1.54 -5.13
## TRAPLINE3 -0.79 1.39 -3.47
## TRAPLINE5 -3.35 1.44 -6.24
## TRAPLINE7 -0.65 1.42 -3.38
## TRAPLINE8 -5.54 1.19 -7.95
## YEAR2 -0.05 0.38 -0.80
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 1.80 0.42 1.01
## SPECIES1:TRAPLINE1 2.27 1.91 -1.35
## SPECIES3:TRAPLINE1 3.28 2.08 -0.73
## SPECIES1:TRAPLINE2 -0.71 1.70 -4.07
## SPECIES3:TRAPLINE2 -0.32 2.02 -4.28
## SPECIES1:TRAPLINE3 -0.70 1.53 -3.80
## SPECIES3:TRAPLINE3 -0.35 2.12 -4.59
## SPECIES1:TRAPLINE5 -0.79 1.66 -4.16
## SPECIES3:TRAPLINE5 -1.56 1.95 -5.47
## SPECIES1:TRAPLINE7 -1.62 1.57 -4.74
## SPECIES3:TRAPLINE7 -2.32 1.97 -6.40
## SPECIES1:TRAPLINE8 1.43 1.44 -1.50
## SPECIES3:TRAPLINE8 3.35 1.68 0.12
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 0.05 0.45 -0.83
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 0.27 0.54 -0.82
## u-95% CI Rhat Bulk_ESS
## Intercept 9.23 1.00 1201
## SPECIES1 -2.65 1.00 1214
## SPECIES3 -1.64 1.00 1450
## TRAPLINE1 -3.00 1.00 1774
## TRAPLINE2 0.84 1.00 1479
## TRAPLINE3 2.00 1.00 1571
## TRAPLINE5 -0.54 1.00 1627
## TRAPLINE7 2.13 1.00 1473
## TRAPLINE8 -3.29 1.00 1335
## YEAR2 0.69 1.00 1817
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 2.66 1.00 1573
## SPECIES1:TRAPLINE1 6.05 1.00 1481
## SPECIES3:TRAPLINE1 7.30 1.00 1644
## SPECIES1:TRAPLINE2 2.49 1.00 1551
## SPECIES3:TRAPLINE2 3.58 1.00 1546
## SPECIES1:TRAPLINE3 2.26 1.00 1312
## SPECIES3:TRAPLINE3 3.70 1.00 1761
## SPECIES1:TRAPLINE5 2.52 1.00 1590
## SPECIES3:TRAPLINE5 2.33 1.00 1713
## SPECIES1:TRAPLINE7 1.43 1.00 1585
## SPECIES3:TRAPLINE7 1.40 1.00 1657
## SPECIES1:TRAPLINE8 4.20 1.00 1507
## SPECIES3:TRAPLINE8 6.62 1.00 1474
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 0.92 1.00 1521
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 1.34 1.00 1681
## Tail_ESS
## Intercept 1620
## SPECIES1 1426
## SPECIES3 1749
## TRAPLINE1 1684
## TRAPLINE2 1617
## TRAPLINE3 1775
## TRAPLINE5 1722
## TRAPLINE7 1787
## TRAPLINE8 1398
## YEAR2 1466
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 1898
## SPECIES1:TRAPLINE1 1538
## SPECIES3:TRAPLINE1 1655
## SPECIES1:TRAPLINE2 1821
## SPECIES3:TRAPLINE2 1718
## SPECIES1:TRAPLINE3 1673
## SPECIES3:TRAPLINE3 1873
## SPECIES1:TRAPLINE5 1505
## SPECIES3:TRAPLINE5 1714
## SPECIES1:TRAPLINE7 1853
## SPECIES3:TRAPLINE7 1682
## SPECIES1:TRAPLINE8 1659
## SPECIES3:TRAPLINE8 1810
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 1844
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 1895
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).joint tests
post-hoc analysis
comm.comp1.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="TRAPLINE") |> summary()comm.comp1.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="SPECIES") |> summary()probabilities
## SPECIES TRAPLINE response lower.HPD upper.HPD
## 2 4 0.9991 0.995720 1.0000
## 1 4 0.8955 0.682984 0.9875
## 3 4 0.9449 0.732562 0.9984
## 2 1 0.6473 0.100637 0.9995
## 1 1 0.1224 0.037927 0.2236
## 3 1 0.4401 0.133342 0.7608
## 2 2 0.9912 0.939314 1.0000
## 1 2 0.3215 0.133179 0.5222
## 3 2 0.5890 0.129636 0.9408
## 2 3 0.9978 0.987916 1.0000
## 1 3 0.6614 0.413431 0.8758
## 3 3 0.8441 0.313672 0.9999
## 2 5 0.9729 0.841554 0.9998
## 1 5 0.1232 0.024561 0.2876
## 3 5 0.1153 0.003966 0.4342
## 2 7 0.9982 0.988212 1.0000
## 1 7 0.4710 0.249584 0.7214
## 3 7 0.4781 0.090867 0.8773
## 2 8 0.8112 0.545423 0.9767
## 1 8 0.1238 0.018941 0.3104
## 3 8 0.6581 0.323746 0.9399
##
## Results are averaged over the levels of: YEAR
## Point estimate displayed: median
## Results are back-transformed from the logit scale
## HPD interval probability: 0.95prob. comparisons
mtsqst1 <- comm.comp1.2025 |> emmeans(pairwise~SPECIES*TRAPLINE)
mtsqrt2 <- mtsqst1$contrasts |> gather_emmeans_draws() |>
filter(contrast %in% c(
"SPECIES2 TRAPLINE4 - SPECIES1 TRAPLINE4",
"SPECIES2 TRAPLINE4 - SPECIES3 TRAPLINE4",
"SPECIES2 TRAPLINE2 - SPECIES1 TRAPLINE2",
"SPECIES2 TRAPLINE2 - SPECIES3 TRAPLINE2",
"SPECIES2 TRAPLINE3 - SPECIES1 TRAPLINE3",
"SPECIES2 TRAPLINE3 - SPECIES3 TRAPLINE3",
"SPECIES2 TRAPLINE5 - SPECIES1 TRAPLINE5",
"SPECIES2 TRAPLINE5 - SPECIES3 TRAPLINE5",
"SPECIES2 TRAPLINE7 - SPECIES1 TRAPLINE7",
"SPECIES2 TRAPLINE7 - SPECIES3 TRAPLINE7",
"SPECIES2 TRAPLINE8 - SPECIES1 TRAPLINE8",
"SPECIES2 TRAPLINE8 - SPECIES3 TRAPLINE8",
"SPECIES2 TRAPLINE1 - SPECIES1 TRAPLINE1",
"SPECIES2 TRAPLINE1 - SPECIES3 TRAPLINE1"
))
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n()) 
summary figure
var <- get_variables(comm.comp1.2025)
vari <- get_variables(comm.comp1.2025)[c(1,4:9,12:23)]
comm.comp1_int_draws2 <- comm.comp1.2025 |> spread_draws(!!!syms(vari))
comm.comp1_draws <- comm.comp1.2025 |> gather_draws(b_Intercept, b_SPECIES1, b_SPECIES3) |>
left_join(comm.comp1_int_draws2, by = c(".chain",".iteration",".draw")) |>
mutate(TRAPLINE4_mean = case_when(`.variable` == 'b_Intercept' ~
`.value`,
`.variable` != 'b_Intercept' ~
`.value` + b_Intercept),
TRAPLINE1_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE1,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE1 + `b_SPECIES1:TRAPLINE1`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE1 + `b_SPECIES3:TRAPLINE1`),
TRAPLINE2_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE2,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE2 + `b_SPECIES1:TRAPLINE2`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE2 + `b_SPECIES3:TRAPLINE2`),
TRAPLINE3_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE3,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE3 + `b_SPECIES1:TRAPLINE3`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE3 + `b_SPECIES3:TRAPLINE3`),
TRAPLINE5_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE5,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE5 + `b_SPECIES1:TRAPLINE5`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE5 + `b_SPECIES3:TRAPLINE5`),
TRAPLINE7_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE7,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE7 + `b_SPECIES1:TRAPLINE7`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE7 + `b_SPECIES3:TRAPLINE7`),
TRAPLINE8_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE8,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE8 + `b_SPECIES1:TRAPLINE8`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE8 + `b_SPECIES3:TRAPLINE8`),
transformed_trapline4 = exp(TRAPLINE4_mean)/(1+exp(TRAPLINE4_mean)),
transformed_trapline1 = exp(TRAPLINE1_mean)/(1+exp(TRAPLINE1_mean)),
transformed_trapline2 = exp(TRAPLINE2_mean)/(1+exp(TRAPLINE2_mean)),
transformed_trapline3 = exp(TRAPLINE3_mean)/(1+exp(TRAPLINE3_mean)),
transformed_trapline5 = exp(TRAPLINE5_mean)/(1+exp(TRAPLINE5_mean)),
transformed_trapline7 = exp(TRAPLINE7_mean)/(1+exp(TRAPLINE7_mean)),
transformed_trapline8 = exp(TRAPLINE8_mean)/(1+exp(TRAPLINE8_mean)))
comm.comp1_draws_plotting <- comm.comp1_draws |>
pivot_longer(cols = starts_with("transformed"),
names_to = "TRAPLINE",
values_to = "infect_prob") %>%
transmute(species = case_when(`.variable` == "b_Intercept" ~ "RBV",
`.variable` == "b_SPECIES1" ~ "DM",
`.variable` == "b_SPECIES3" ~ "WJM"),
trapline = case_when(TRAPLINE == 'transformed_trapline4' ~ "trapline4",
TRAPLINE == 'transformed_trapline1' ~ "trapline1",
TRAPLINE == 'transformed_trapline2' ~ "trapline2",
TRAPLINE == 'transformed_trapline3' ~ "trapline3",
TRAPLINE == 'transformed_trapline5' ~ "trapline5",
TRAPLINE == 'transformed_trapline7' ~ "trapline7",
TRAPLINE == 'transformed_trapline8' ~ "trapline8"),
infect_prob = infect_prob,
speciesb = species,
traplineb = trapline,
chain = `.chain`,
iteration = `.iteration`,
draw_num = `.draw`) %>%
unite("id",speciesb,traplineb,sep = "_") %>%
mutate(id = factor(id, levels = c("RBV_trapline4","DM_trapline4","WJM_trapline4",
"RBV_trapline1","DM_trapline1","WJM_trapline1",
"RBV_trapline2","DM_trapline2","WJM_trapline2",
"RBV_trapline3","DM_trapline3","WJM_trapline3",
"RBV_trapline5","DM_trapline5","WJM_trapline5",
"RBV_trapline7","DM_trapline7","WJM_trapline7",
"RBV_trapline8","DM_trapline8","WJM_trapline8")))
trapline_names <- c("Sugar maple hardwood", " Cut-over mixed-wood", "Dense mixed-wood",
"Conifer", "White pine/white spruce", "Black spruce/aspen",
"White/red pine")
trapline_index <- c("trapline1","trapline2","trapline3","trapline4","trapline5",
"trapline7","trapline8")
species_labels <- c("DM" = "Deer mouse","RBV" = "Red-Backed Vole","WJM" = "Woodland Jumping Mouse") ridge_plot1 <- comm.comp1_draws_plotting |>
ggplot(aes(infect_prob, trapline, fill = species)) +
scale_fill_manual(values = alpha(c("#7C8081","#7E4640", "#ABA159"),0.7)) +
geom_density_ridges_gradient(scale=2, rel_min_height = 0.01) +
facet_wrap(~species, labeller = labeller(species = species_labels))+
theme_bw()+
scale_x_continuous(breaks = seq(0, 1, by = .2))+
xlab("Probability of infection")+ylab("Habitat type")+
scale_y_discrete(limits =c("trapline1",
"trapline8",
"trapline5",
"trapline2",
"trapline7",
"trapline3",
"trapline4"),
labels = c("Sugar maple hardwood",
"White/red pine",
"White pine/white spruce",
"Cut-over mixed-wood",
"Black spruce/aspen",
"Dense mixed-wood",
"Conifer"))+
theme(axis.text=element_text(size=18),
axis.title=element_text(size=25,face="bold"),
strip.text = element_text(size=18),
legend.position = "none",
panel.grid.major = element_blank(), panel.grid.minor = element_blank());ridge_plot1## Picking joint bandwidth of 0.018## Picking joint bandwidth of 0.0129## Picking joint bandwidth of 0.0332
summary figure 2
date.plot <- orangemitedata_y1y2 |>
group_by(SPECIES, TRAPLINE) |>
modelr::data_grid(DATE_JULIAN=seq(from=120,
to=max(orangemitedata_y1y2$DATE_JULIAN),
length.out=100),
YEAR=1) |>
add_fitted_draws(comm.comp1.2025,
n =100,
re_formula=NA) |>
#unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=DATE_JULIAN, color=SPECIES)) +
geom_line(aes(y=.value, group=paste(SPECIES, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=SPECIES), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~TRAPLINE); date.plot## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
## named `object` and the `n` parameter is now named `ndraws`.## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Low-RBV year
Fit models
priors
model formula
Fit model
comm.comp2.2025 <- brm(comm.comp2.formula,
data =orangemitedata_y3y4,
family=bernoulli(),
prior = priors_y3y4,
iter =8000,
warmup = 4000,
seed=123,
thin=2,
chains = 4,
cores = 2,
save_pars = save_pars(all = TRUE),
sample_prior = "yes",
control = list(adapt_delta=0.99, max_treedepth=15))
saveRDS(comm.comp2.2025, "models/comm_comp2_2025.RDS")model validations (convergence)
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model investigation
summary
## Family: bernoulli
## Links: mu = logit
## Formula: ORANGE_MITES ~ SPECIES + TRAPLINE + YEAR + scale(DATE_JULIAN, center = TRUE, scale = TRUE) + SPECIES:TRAPLINE + SPECIES:scale(DATE_JULIAN, center = TRUE, scale = TRUE) + (1 | TAG_LEFT)
## Data: orangemitedata_y3y4 (Number of observations: 1125)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 2;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~TAG_LEFT (Number of levels: 436)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 4.46 0.69 3.27 5.96 1.00 4538 6076
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI
## Intercept 4.97 1.73 1.65
## SPECIES1 -8.29 1.89 -12.13
## SPECIES3 -11.03 3.22 -17.68
## TRAPLINE1 -3.43 3.69 -10.71
## TRAPLINE2 -0.34 2.40 -5.01
## TRAPLINE3 -1.31 2.70 -6.58
## TRAPLINE5 -1.11 2.13 -5.19
## TRAPLINE7 -0.12 2.07 -4.25
## TRAPLINE8 -3.95 3.06 -10.02
## YEAR4 -2.84 0.84 -4.63
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 1.08 0.59 -0.03
## SPECIES1:TRAPLINE1 -3.40 3.72 -10.61
## SPECIES3:TRAPLINE1 0.00 4.96 -9.85
## SPECIES1:TRAPLINE2 -0.26 2.51 -5.15
## SPECIES3:TRAPLINE2 -0.09 5.04 -9.87
## SPECIES1:TRAPLINE3 -1.58 2.84 -7.19
## SPECIES3:TRAPLINE3 0.04 4.99 -9.65
## SPECIES1:TRAPLINE5 -1.26 2.40 -6.03
## SPECIES3:TRAPLINE5 -2.23 4.29 -10.98
## SPECIES1:TRAPLINE7 -1.98 2.27 -6.61
## SPECIES3:TRAPLINE7 -4.47 3.77 -11.96
## SPECIES1:TRAPLINE8 1.10 3.18 -5.05
## SPECIES3:TRAPLINE8 -2.27 4.21 -10.77
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE -0.50 0.63 -1.75
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE -1.21 2.47 -6.11
## u-95% CI Rhat Bulk_ESS
## Intercept 8.48 1.00 6528
## SPECIES1 -4.72 1.00 6218
## SPECIES3 -5.03 1.00 6833
## TRAPLINE1 3.77 1.00 7464
## TRAPLINE2 4.27 1.00 6261
## TRAPLINE3 4.00 1.00 7257
## TRAPLINE5 3.02 1.00 5663
## TRAPLINE7 4.00 1.00 6312
## TRAPLINE8 1.79 1.00 7430
## YEAR4 -1.34 1.00 6184
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 2.25 1.00 7013
## SPECIES1:TRAPLINE1 3.84 1.00 7646
## SPECIES3:TRAPLINE1 9.64 1.00 6523
## SPECIES1:TRAPLINE2 4.74 1.00 6063
## SPECIES3:TRAPLINE2 9.82 1.00 6092
## SPECIES1:TRAPLINE3 3.96 1.00 6763
## SPECIES3:TRAPLINE3 9.44 1.00 6022
## SPECIES1:TRAPLINE5 3.35 1.00 5956
## SPECIES3:TRAPLINE5 5.90 1.00 6866
## SPECIES1:TRAPLINE7 2.42 1.00 6431
## SPECIES3:TRAPLINE7 2.69 1.00 7081
## SPECIES1:TRAPLINE8 7.32 1.00 7088
## SPECIES3:TRAPLINE8 5.70 1.00 6517
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 0.69 1.00 6881
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 3.70 1.00 7080
## Tail_ESS
## Intercept 7175
## SPECIES1 6756
## SPECIES3 7320
## TRAPLINE1 6592
## TRAPLINE2 6549
## TRAPLINE3 7339
## TRAPLINE5 6458
## TRAPLINE7 7181
## TRAPLINE8 6645
## YEAR4 6984
## scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 7252
## SPECIES1:TRAPLINE1 6686
## SPECIES3:TRAPLINE1 6019
## SPECIES1:TRAPLINE2 6949
## SPECIES3:TRAPLINE2 5818
## SPECIES1:TRAPLINE3 6885
## SPECIES3:TRAPLINE3 6171
## SPECIES1:TRAPLINE5 6406
## SPECIES3:TRAPLINE5 6962
## SPECIES1:TRAPLINE7 6738
## SPECIES3:TRAPLINE7 7366
## SPECIES1:TRAPLINE8 6872
## SPECIES3:TRAPLINE8 6773
## SPECIES1:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 7250
## SPECIES3:scaleDATE_JULIANcenterEQTRUEscaleEQTRUE 7422
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).joint tests
post-hoc analysis
comm.comp2.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="TRAPLINE") |> summary()comm.comp2.2025 |> emmeans(~ SPECIES | TRAPLINE, type ="response") |> regrid() |> pairs(by="SPECIES") |> summary()probabilities
## SPECIES TRAPLINE response lower.HPD upper.HPD
## 2 4 9.71e-01 7.03e-01 0.999963
## 1 4 9.05e-03 1.71e-05 0.056271
## 3 4 6.56e-04 0.00e+00 0.099588
## 2 1 5.31e-01 1.50e-03 0.999997
## 1 1 1.14e-05 0.00e+00 0.000206
## 3 1 1.83e-05 0.00e+00 0.660663
## 2 2 9.60e-01 4.14e-01 0.999997
## 1 2 5.19e-03 2.36e-05 0.023954
## 3 2 3.82e-04 0.00e+00 0.920104
## 2 3 8.99e-01 1.14e-01 0.999998
## 1 3 5.53e-04 4.00e-07 0.004829
## 3 3 1.61e-04 0.00e+00 0.858584
## 2 5 9.17e-01 4.00e-01 0.999993
## 1 5 9.24e-04 3.00e-07 0.007364
## 3 5 2.71e-05 0.00e+00 0.031493
## 2 7 9.66e-01 6.97e-01 0.999951
## 1 7 1.16e-03 2.00e-06 0.007082
## 3 7 7.60e-06 0.00e+00 0.001812
## 2 8 4.15e-01 4.70e-06 0.992098
## 1 8 5.67e-04 0.00e+00 0.005284
## 3 8 1.80e-06 0.00e+00 0.001603
##
## Results are averaged over the levels of: YEAR
## Point estimate displayed: median
## Results are back-transformed from the logit scale
## HPD interval probability: 0.95mtsqst1 <- comm.comp2.2025 |> emmeans(pairwise~SPECIES*TRAPLINE)
mtsqrt2 <- mtsqst1$contrasts |> gather_emmeans_draws() |>
filter(contrast %in% c(
"SPECIES2 TRAPLINE4 - SPECIES1 TRAPLINE4",
#"SPECIES2 TRAPLINE4 - SPECIES3 TRAPLINE4",
"SPECIES2 TRAPLINE2 - SPECIES1 TRAPLINE2",
#"SPECIES2 TRAPLINE2 - SPECIES3 TRAPLINE2",
"SPECIES2 TRAPLINE3 - SPECIES1 TRAPLINE3",
#"SPECIES2 TRAPLINE3 - SPECIES3 TRAPLINE3",
#"SPECIES2 TRAPLINE5 - SPECIES1 TRAPLINE5",
#"SPECIES2 TRAPLINE5 - SPECIES3 TRAPLINE5",
"SPECIES2 TRAPLINE7 - SPECIES1 TRAPLINE7",
#"SPECIES2 TRAPLINE7 - SPECIES3 TRAPLINE7",
"SPECIES2 TRAPLINE8 - SPECIES1 TRAPLINE8"
#"SPECIES2 TRAPLINE8 - SPECIES3 TRAPLINE8",
#"SPECIES2 TRAPLINE1 - SPECIES1 TRAPLINE1",
#"SPECIES2 TRAPLINE1 - SPECIES3 TRAPLINE1"
))
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n()) 
summary figure
var <- get_variables(comm.comp2.2025)
vari <- get_variables(comm.comp2.2025)[c(1,4:9,12:23)]
comm.comp2_int_draws2 <- comm.comp2.2025 |> spread_draws(!!!syms(vari))
comm.comp2_draws <- comm.comp2.2025 |> gather_draws(b_Intercept, b_SPECIES1, b_SPECIES3) |>
left_join(comm.comp2_int_draws2, by = c(".chain",".iteration",".draw")) |>
mutate(TRAPLINE4_mean = case_when(`.variable` == 'b_Intercept' ~
`.value`,
`.variable` != 'b_Intercept' ~
`.value` + b_Intercept),
TRAPLINE1_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE1,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE1 + `b_SPECIES1:TRAPLINE1`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE1 + `b_SPECIES3:TRAPLINE1`),
TRAPLINE2_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE2,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE2 + `b_SPECIES1:TRAPLINE2`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE2 + `b_SPECIES3:TRAPLINE2`),
TRAPLINE3_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE3,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE3 + `b_SPECIES1:TRAPLINE3`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE3 + `b_SPECIES3:TRAPLINE3`),
TRAPLINE5_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE5,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE5 + `b_SPECIES1:TRAPLINE5`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE5 + `b_SPECIES3:TRAPLINE5`),
TRAPLINE7_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE7,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE7 + `b_SPECIES1:TRAPLINE7`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE7 + `b_SPECIES3:TRAPLINE7`),
TRAPLINE8_mean = case_when(`.variable` == 'b_Intercept' ~
`.value` + b_TRAPLINE8,
`.variable` == 'b_SPECIES1' ~
`.value` + b_Intercept +
b_TRAPLINE8 + `b_SPECIES1:TRAPLINE8`,
`.variable` == 'b_SPECIES3' ~
`.value` + b_Intercept +
b_TRAPLINE8 + `b_SPECIES3:TRAPLINE8`),
transformed_trapline4 = exp(TRAPLINE4_mean)/(1+exp(TRAPLINE4_mean)),
transformed_trapline1 = exp(TRAPLINE1_mean)/(1+exp(TRAPLINE1_mean)),
transformed_trapline2 = exp(TRAPLINE2_mean)/(1+exp(TRAPLINE2_mean)),
transformed_trapline3 = exp(TRAPLINE3_mean)/(1+exp(TRAPLINE3_mean)),
transformed_trapline5 = exp(TRAPLINE5_mean)/(1+exp(TRAPLINE5_mean)),
transformed_trapline7 = exp(TRAPLINE7_mean)/(1+exp(TRAPLINE7_mean)),
transformed_trapline8 = exp(TRAPLINE8_mean)/(1+exp(TRAPLINE8_mean)))
comm.comp2_draws_plotting <- comm.comp2_draws %>%
pivot_longer(cols = starts_with("transformed"),
names_to = "TRAPLINE",
values_to = "infect_prob") %>%
transmute(species = case_when(`.variable` == "b_Intercept" ~ "RBV",
`.variable` == "b_SPECIES1" ~ "DM",
`.variable` == "b_SPECIES3" ~ "WJM"),
trapline = case_when(TRAPLINE == 'transformed_trapline4' ~ "trapline4",
TRAPLINE == 'transformed_trapline1' ~ "trapline1",
TRAPLINE == 'transformed_trapline2' ~ "trapline2",
TRAPLINE == 'transformed_trapline3' ~ "trapline3",
TRAPLINE == 'transformed_trapline5' ~ "trapline5",
TRAPLINE == 'transformed_trapline7' ~ "trapline7",
TRAPLINE == 'transformed_trapline8' ~ "trapline8"),
infect_prob = infect_prob,
speciesb = species,
traplineb = trapline,
chain = `.chain`,
iteration = `.iteration`,
draw_num = `.draw`) %>%
unite("id",speciesb,traplineb,sep = "_") %>%
mutate(id = factor(id, levels = c("RBV_trapline4","DM_trapline4","WJM_trapline4",
"RBV_trapline1","DM_trapline1","WJM_trapline1",
"RBV_trapline2","DM_trapline2","WJM_trapline2",
"RBV_trapline3","DM_trapline3","WJM_trapline3",
"RBV_trapline5","DM_trapline5","WJM_trapline5",
"RBV_trapline7","DM_trapline7","WJM_trapline7",
"RBV_trapline8","DM_trapline8","WJM_trapline8")))
trapline_names <- c("Sugar maple hardwood", " Cut-over mixed-wood", "Dense mixed-wood",
"Conifer", "White pine/white spruce", "Black spruce/aspen",
"White/red pine")
trapline_index <- c("trapline1","trapline2","trapline3","trapline4","trapline5",
"trapline7","trapline8")
species_labels <- c("DM" = "Deer mouse","RBV" = "Red-Backed Vole","WJM" = "Woodland Jumping Mouse") ridge_plot2 <- comm.comp2_draws_plotting %>%
ggplot(aes(infect_prob, trapline, fill = species)) +
scale_fill_manual(values = alpha(c("#7C8081","#7E4640", "#ABA159"),0.7)) +
geom_density_ridges_gradient(scale = 2, rel_min_height = 0.01)+
facet_wrap(~species, labeller = labeller(species = species_labels))+
theme_bw()+
scale_x_continuous(breaks = seq(0, 1, by = .2))+
xlab("P(Infection)")+ylab("Habitat type")+
scale_y_discrete(limits =c("trapline1",
"trapline8",
"trapline5",
"trapline2",
"trapline7",
"trapline3",
"trapline4"),
labels = c("Sugar maple hardwood",
"White/red pine",
"White pine/white spruce",
"Cut-over mixed-wood",
"Black spruce/aspen",
"Dense mixed-wood",
"Conifer"))+
theme(axis.text=element_text(size=18),
axis.title=element_text(size=25,face="bold"),
strip.text = element_text(size=18),
legend.position = "none",
panel.grid.major = element_blank(), panel.grid.minor = element_blank());ridge_plot2## Picking joint bandwidth of 0.00185## Picking joint bandwidth of 0.0204## Picking joint bandwidth of 0.00299
summary figure 2
date.plot <- orangemitedata_y3y4 |>
group_by(SPECIES, TRAPLINE) |>
modelr::data_grid(DATE_JULIAN=seq(from=120,
to=max(orangemitedata_y3y4$DATE_JULIAN),
length.out=100),
YEAR=3) |>
add_fitted_draws(comm.comp2.2025,
n =100,
re_formula=NA) |>
#unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=DATE_JULIAN, color=SPECIES)) +
geom_line(aes(y=.value, group=paste(SPECIES, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=SPECIES), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~TRAPLINE); date.plot## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
## means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
## named `object` and the `n` parameter is now named `ndraws`.## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
RBV community presence
Fit models
Data manipulation
priors
model formula
Fit model
rbv.perc.2025 <- brm(rbv.perc_formula,
data =orangemitedata_perc.RBV,
family=bernoulli(),
prior =rbv.perc_priors,
iter =2000,
warmup =1000,
seed=123,
thin=2,
chains = 4,
cores = 2,
save_pars = save_pars(all = TRUE),
sample_prior = "yes",
control = list(adapt_delta=0.99, max_treedepth=15))
saveRDS(rbv.perc.2025, "models/rbv_perc_2025.RDS")model validations (convergence)
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model investigation
summary
## Family: bernoulli
## Links: mu = logit
## Formula: ORANGE_MITES ~ PERC_RBV * SPECIES + (1 | TAG_LEFT)
## Data: orangemitedata_perc.RBV (Number of observations: 2571)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~TAG_LEFT (Number of levels: 1037)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.01 0.24 2.57 3.53 1.01 907 1373
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.51 0.61 0.38 2.70 1.00 1190 1629
## PERC_RBV 8.24 1.72 4.98 11.66 1.00 1411 1606
## SPECIES1 -5.02 0.69 -6.44 -3.71 1.01 1068 1407
## SPECIES3 -3.85 0.78 -5.35 -2.38 1.01 1382 1828
## PERC_RBV:SPECIES1 -2.81 1.95 -6.64 0.84 1.00 1401 1441
## PERC_RBV:SPECIES3 1.33 2.33 -3.24 6.14 1.00 1478 1753
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).hpd.summary
post-hoc analysis
rbv.perc.2025 |> emtrends(var = "PERC_RBV", type = "response") |> regrid() |> summary(by = NULL, adjust = "tukey", infer=TRUE) rbv.perc.2025 |> emtrends(var = "PERC_RBV", type = "response") |>
pairs(by = "PERC_RBV") |>
summary(by = NULL, adjust = "tukey", infer=TRUE) ## SPECIES PERC_RBV response lower.HPD upper.HPD
## 2 0.156 0.9412 0.8866 0.9810
## 1 0.156 0.0652 0.0402 0.0961
## 3 0.156 0.3027 0.1647 0.4472
##
## Point estimate displayed: median
## Results are back-transformed from the logit scale
## HPD interval probability: 0.95rbv.perc.2025 |> emmeans(~SPECIES*PERC_RBV, type = "response") |>
pairs(by = "PERC_RBV") |>
summary(by = NULL,
adjust = "tukey",
infer=TRUE) 
Summary
perc_rbv.figure <- orangemitedata_v3 |>
group_by(SPECIES) |>
modelr::data_grid(PERC_RBV=seq(from = min(orangemitedata_v3$PERC_RBV),
to=max(orangemitedata_v3$PERC_RBV),
length.out = 101)) |>
add_linpred_draws(rbv.perc.2025, ndraws=100,
re_formula = NA) |>
mutate(.linpred =exp(.linpred)/(1+exp(.linpred))) |>
ggplot(aes(x=PERC_RBV, color=SPECIES))+
geom_line(aes(y=.linpred, group = paste(SPECIES, .draw)), alpha=0.2)+
geom_smooth(aes(y=.linpred, color=SPECIES), se=F, size=1.5)+theme_classic()+
xlab("% red backed vole in community composition")+ylab("Probability of infection")+
scale_color_manual(labels=c("Red backed vole", "Deer mouse", "Woodland jumping mouse"),
values=c("#7E4640","#7C8081","#ABA159")) +
theme(legend.position = "bottom",
legend.title = element_blank(),
axis.title.x = element_text(size=25),
axis.text.x = element_text(size=18),
axis.title.y = element_text(size=25),
axis.text.y = element_text(size=18),
legend.text = element_text(size=15),
axis.line = element_line(size=1.2)) ## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Body condition
Data manipulation
## Rows: 51 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (2): julian_date, sampling_period
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.bc <- read_excel("data/Small Mammal Project 2016 (ALL)_cleaned.xlsx",
sheet = "Falls Lines") |>
clean_names() |> #use the janitor package to clean up the names so everything's consistent
filter(age == 'A') |> #subset to only adults JV
filter(!(repro_cond == 'PERF/PREG' | repro_cond == 'Preg' | repro_cond == 'PREG' |
repro_cond == 'Preg/Lac'| repro_cond == 'PREG/LAC' |
repro_cond == 'PREG/PERF')) |> #subset to remove preg JV
mutate(skull_length = as.numeric(skull_length), tag_left = as.factor(tag_left)) |> #convert skull length to numeric like the rest of the variables JV
filter(species == 'DM' | species == 'RBV') |> #subset the df to only the species we're interested in
mutate(julian_date = yday(date)) |> #group into sampling periods based on date
#connect the sampling csv with this df by matching the julian dates
#from each together
left_join(sampling, by = "julian_date") |>
#so now we have a whole column of sampling periods based on the julian date
#for each observation
#now we order by individual
arrange(tag_left) |>
#the group by individual and sampling period
group_by(tag_left, sampling_period) |>
##create average for weight JV
mutate(avg_weight_sampling_period = mean(weight, na.rm = TRUE)) |>
mutate(tag_left_num = as.numeric(tag_left)) |>
distinct(avg_weight_sampling_period, tag_left_num, .keep_all = TRUE) |>
#now remove that grouping variable so we can group by individual over the
#whole summer now
ungroup() |>
#by individual across all samples
group_by(tag_left) |>
#make a new column to average measures we care about JV
mutate(avg_tail_summer = mean(tail, na.rm = TRUE),
avg_hindfoot_summer = mean(hindfoot, na.rm = TRUE),
avg_ear_summer = mean(ear, na.rm = TRUE),
avg_skull_width_summer = mean(skull_width, na.rm = TRUE),
avg_skull_length_summer = mean(skull_length, na.rm = TRUE)) |>
##creating infection stages col JV
mutate(group = case_when(any(orange_mites_y_n == "Y") ~ "Parasitized", TRUE ~ "Uninfected")) |>
ungroup() |>
mutate(infection_stage = case_when((group == "Parasitized" & orange_mites_y_n == "Y") ~ "After",
(group == "Parasitized" & orange_mites_y_n == "N") ~ "Before",
(group == "Uninfected") ~ "Before",
TRUE ~ "NA")) |>
##There were two individuals that switched from having mites to not having mites
##The code below makes sure they are still listed as "after" when they lose their mites JV
mutate(infection_stage = case_when((tag_left == 50468 & julian_date > 196) ~ "After",
(tag_left == 50968 & julian_date > 207) ~ "After",
TRUE ~ infection_stage)) |>
mutate(group = as.factor(group), infection_stage = as.factor(infection_stage)) |>
filter(infection_stage != "NA") |>
unite(infection, c("group", "infection_stage"), remove = FALSE) |>
mutate(infection = as.factor(infection)) |>
##removing NA's JV
drop_na(avg_skull_width_summer, avg_skull_length_summer, avg_hindfoot_summer,
avg_weight_sampling_period) |>
##formatting predictors JV
mutate(sex = as.factor(sex)) Split dataframes
Red-backed voles
Princple component analysis (PCA)
pca_rbv <- princomp(~ avg_hindfoot_summer + avg_skull_length_summer + avg_skull_width_summer,
rbv,
cor = TRUE)
loadings(pca_rbv)##
## Loadings:
## Comp.1 Comp.2 Comp.3
## avg_hindfoot_summer 0.562 0.747 0.355
## avg_skull_length_summer 0.599 -0.798
## avg_skull_width_summer 0.570 -0.661 0.487
##
## Comp.1 Comp.2 Comp.3
## SS loadings 1.000 1.000 1.000
## Proportion Var 0.333 0.333 0.333
## Cumulative Var 0.333 0.667 1.000## Importance of components:
## Comp.1 Comp.2 Comp.3
## Standard deviation 1.3071989 0.8296390 0.7764858
## Proportion of Variance 0.5695896 0.2294336 0.2009768
## Cumulative Proportion 0.5695896 0.7990232 1.0000000
## Comp.1 Comp.2 Comp.3
## 1.7087689 0.6883008 0.6029303
Extracting residuals
bodycond_correlation_rbv <- lm(avg_weight_sampling_period ~ body_size_rbv,
data =rbv)
rbv$bodycondition <- bodycond_correlation_rbv$residuals
hist(rbv$bodycondition)
priors
model formula
Fit model (work in progress)
RBVbodycondition <- brm(RBVbodycondition_formula,
data =rbv,
family=bernoulli(),
prior =RBVbodycondition_priors,
iter =2000,
warmup =1000,
seed=123,
thin=2,
chains = 4,
cores = 2,
save_pars = save_pars(all = TRUE),
sample_prior = "yes",
control = list(adapt_delta=0.99, max_treedepth=15))
vsaveRDS(RBVbodycondition, "models/RBVbodycondition_2025.RDS")model validations (convergence)
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model investigation
summary
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: bodycondition ~ infection + sex * scale(julian_date, center = TRUE, scale = TRUE) + (1 | tag_left)
## Data: rbv (Number of observations: 115)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~tag_left (Number of levels: 86)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.18 0.37 2.49 3.93 1.01 588 1022
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI
## Intercept 1.08 0.68 -0.21
## infectionParasitized_Before -0.96 1.13 -3.14
## infectionUninfected_Before 0.44 1.18 -1.88
## sexM -2.05 0.83 -3.63
## scalejulian_datecenterEQTRUEscaleEQTRUE 0.98 0.60 -0.25
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE -0.95 0.69 -2.31
## u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.45 1.00 1074 1506
## infectionParasitized_Before 1.25 1.00 1581 1743
## infectionUninfected_Before 2.68 1.00 1415 1650
## sexM -0.44 1.00 1195 1465
## scalejulian_datecenterEQTRUEscaleEQTRUE 2.10 1.00 1434 1617
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE 0.37 1.00 1584 1606
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.88 0.28 1.42 2.49 1.01 553 852
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).posteriors <- as.data.frame(RBVbodycondition)
covariate_samples <- posteriors |> select(b_Intercept, b_infectionParasitized_Before,b_infectionUninfected_Before)
combined_effect <- rowSums(covariate_samples)
hdi(combined_effect, credMass =0.95)## lower upper
## -2.958808 4.128201
## attr(,"credMass")
## [1] 0.95post-hoc analysis
Infection
## infection emmean lower.HPD upper.HPD
## Parasitized_After 0.0521 -0.827 0.962
## Parasitized_Before -0.8898 -3.117 1.222
## Uninfected_Before 0.5414 -1.525 2.531
##
## Results are averaged over the levels of: sex
## Point estimate displayed: median
## HPD interval probability: 0.95RBVbodycondition |> emmeans(~infection, type = "response") |>
pairs() |> summary(by = NULL, adjust = "tukey", infer=TRUE) sex
## NOTE: Results may be misleading due to involvement in interactions## sex emmean lower.HPD upper.HPD
## F 0.921 -0.708 2.482
## M -1.144 -2.311 0.147
##
## Results are averaged over the levels of: infection
## Point estimate displayed: median
## HPD interval probability: 0.95RBVbodycondition |> emmeans(~sex, type = "response") |>
pairs() |>
summary(by = NULL,
adjust = "tukey",
infer=TRUE) ## NOTE: Results may be misleading due to involvement in interactions
summary figure
Deer mice
Princple component analysis (PCA)
pca_dm <- princomp(~ avg_hindfoot_summer + avg_skull_length_summer + avg_skull_width_summer,
dm,
cor = TRUE)
loadings(pca_dm)##
## Loadings:
## Comp.1 Comp.2 Comp.3
## avg_hindfoot_summer 0.656 0.748
## avg_skull_length_summer 0.581 0.567 -0.585
## avg_skull_width_summer 0.482 -0.818 -0.314
##
## Comp.1 Comp.2 Comp.3
## SS loadings 1.000 1.000 1.000
## Proportion Var 0.333 0.333 0.333
## Cumulative Var 0.333 0.667 1.000## Importance of components:
## Comp.1 Comp.2 Comp.3
## Standard deviation 1.1650734 0.9573210 0.8521387
## Proportion of Variance 0.4524654 0.3054879 0.2420468
## Cumulative Proportion 0.4524654 0.7579532 1.0000000
## Comp.1 Comp.2 Comp.3
## 1.3573961 0.9164636 0.7261403
Extracting residuals
bodycond_correlation_dm <- lm(dm$avg_weight_sampling_period ~ body_size_dm)
dm$bodycondition <- bodycond_correlation_dm$residuals
hist(dm$bodycondition)
priors
model formula
Fit model (work in progress)
model validations (convergence)
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
model investigation
summary
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: bodycondition ~ infection + sex * scale(julian_date, center = TRUE, scale = TRUE) + (1 | tag_left)
## Data: dm (Number of observations: 184)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 2;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~tag_left (Number of levels: 88)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.97 0.18 0.63 1.32 1.00 1063 1608
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI
## Intercept 0.06 0.35 -0.63
## infectionParasitized_Before -0.17 0.34 -0.83
## infectionUninfected_Before 0.07 0.35 -0.63
## sexM -0.02 0.32 -0.66
## scalejulian_datecenterEQTRUEscaleEQTRUE -0.71 0.21 -1.11
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE 0.11 0.25 -0.37
## u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.74 1.00 1608 1644
## infectionParasitized_Before 0.48 1.00 1747 1929
## infectionUninfected_Before 0.76 1.00 1519 1726
## sexM 0.60 1.00 1861 1658
## scalejulian_datecenterEQTRUEscaleEQTRUE -0.30 1.00 1877 1788
## sexM:scalejulian_datecenterEQTRUEscaleEQTRUE 0.59 1.00 1817 1771
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.30 0.09 1.13 1.50 1.00 1354 1607
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).post-hoc analysis
Infection
posteriors <- as.data.frame(DMbodycondition)
covariate_samples <- posteriors |> select(b_Intercept, b_infectionParasitized_Before,b_infectionUninfected_Before)
combined_effect <- rowSums(covariate_samples)
hdi(combined_effect, credMass =0.95)## lower upper
## -0.9724559 0.7825951
## attr(,"credMass")
## [1] 0.95## infection emmean lower.HPD upper.HPD
## Parasitized_After 0.0441 -0.513 0.671
## Parasitized_Before -0.1066 -0.737 0.444
## Uninfected_Before 0.1209 -0.263 0.506
##
## Results are averaged over the levels of: sex
## Point estimate displayed: median
## HPD interval probability: 0.95DMbodycondition |> emmeans(~infection, type = "response") |>
pairs() |> summary(by = NULL, adjust = "tukey", infer=TRUE) sex
## NOTE: Results may be misleading due to involvement in interactions## sex emmean lower.HPD upper.HPD
## F 0.03317 -0.467 0.545
## M 0.00735 -0.407 0.426
##
## Results are averaged over the levels of: infection
## Point estimate displayed: median
## HPD interval probability: 0.95DMbodycondition |> emmeans(~sex, type = "response") |>
pairs() |>
summary(by = NULL,
adjust = "tukey",
infer=TRUE) ## NOTE: Results may be misleading due to involvement in interactions
summary figure
Red backed voles
rbv.infect <- as.data.frame(emmeans(RBVbodycondition, ~infection))
rbv.infect.obs <- drop_na(rbv, infection) |>
mutate(Pred =predict(RBVbodycondition, re.form =NULL, type ='response'),
Resid =residuals(RBVbodycondition, type ="ordinary"),
Fit =Pred+Resid) rbv.bodycondition.plot <- ggplot(data =rbv.infect, aes(x=infection, y=emmean)) +
geom_jitter(data =rbv.infect.obs, aes(x= infection, y =Fit[,1]),
width=0.02,
color ="#7E4640",
alpha =0.15) +
geom_pointrange(aes(ymin =lower.HPD,
ymax =upper.HPD),
color ="black")+
xlab("") + ylab("Body condition") +
scale_y_continuous(limits =c(-10,10)) +
scale_x_discrete(limits =c("Uninfected_Before","Parasitized_Before","Parasitized_After")) +
theme_classic() +
theme(legend.position = "none",
axis.text.x = element_blank()); rbv.bodycondition.plot
Deer mice
dm.infect <- as.data.frame(emmeans(DMbodycondition, ~infection))
dm.infect.obs <- drop_na(dm, infection) |>
mutate(Pred =predict(DMbodycondition, re.form =NULL, type ='response'),
Resid =residuals(DMbodycondition, type ="ordinary"),
Fit =Pred+Resid) dm.bodycondition.plot <- ggplot(data =rbv.infect, aes(x=infection, y=emmean)) +
geom_jitter(data =rbv.infect.obs, aes(x= infection, y =Fit[,1]),
width=0.02,
color ="#7C8081",
alpha =0.25) +
geom_pointrange(aes(ymin =lower.HPD,
ymax =upper.HPD),
color ="black") +
xlab("Infection status") + ylab("Body condition") +
scale_y_continuous(limits =c(-10,10)) +
scale_x_discrete(limits =c("Uninfected_Before","Parasitized_Before","Parasitized_After")) +
theme_classic() +
theme(legend.position = "none"); dm.bodycondition.plot
